function [moments,modelSolution,modelSim,xGrid,xDrift,xVolatility] = computeMomentsFromParams(params,algParams)
%COMPUTEMOMENTSFROMPARAMS Computes model moments from set of parameters
%   Inputs:
%     - params: struct array containing all model parameters
%     - algParams: struct array containing all algorithm parameters
%   Outputs:
%     - moments: struct array containing computed model moments
%     - modelSolution: struct array containing model solution functions (on
%                      the state grid)
%     - modelSim: struct array containing simulated paths
%     - xGrid: grid for state variable x:=\tilde{\sigma}^2
%     - xDrift: (arithmetic) drift of state variable on grid
%     - xVolatility: (arithmetic) volatility of state variable on grid
  
  % retrieve parameters
  [sigma,stateVol,stateMeanReversion] = ...
    aux.unpack(params,'sigma','stateVol','stateMeanReversion');
  % retrieve algorithm parameters
  nGridpoints = aux.unpack(algParams,'nGridpoints');

  rng(2387020637); %set random number generator seed

  % (1) construct grids
  stateMean = sigma^2;
  gammaParamA = 2*stateMeanReversion*stateMean/stateVol^2;
  gammaParamB = stateVol^2/2/stateMeanReversion;
  xQuantiles = icdf('Gamma',[0.0005,0.9995],gammaParamA,gammaParamB);
  xGrid = linspace(xQuantiles(1),xQuantiles(2),nGridpoints)';
  sigmaGrid = sqrt(xGrid);
  % x evolution
  xDrift = -stateMeanReversion*(xGrid - stateMean);
  xVolatility = stateVol*sqrt(xGrid);

  % (2) compute the model solution
  modelSolution = solveModelPde(xGrid,sigmaGrid,xDrift,xVolatility,params,algParams);

  % (3) simulate model and compute moments
  [modelSim,moments] = simulateModel(xGrid,xDrift,xVolatility,modelSolution,algParams,true);
               
end

